Boundary Value Problem – Axial Temperature Variation of a current-carrying bare wire¶
(Problem 11.28 from the Numerical Methods for Engineers and Scientists)
This study was developed within the Master’s Degree in Advanced Computing, under the curricular unit Numerical Simulation for Engineering, by Diogo Silva (PG61444) and Tomás Pereira (PG59810).
The objective is to numerically solve the axial temperature distribution along a current-carrying bare wire.
Figure 1 – Axial Temperature Variation of a Current-Carrying Bare Wire
The differential equation that describes the problem is the following:
$$\frac{d^2T}{dx^2} - \frac{4h}{kD}\,(T - T_\infty) - \frac{4 \varepsilon \sigma_{SB}}{kD}\,(T^4 - T_\infty^4) = -\frac{I^2 \rho_e}{k\,(\pi D^2 /4)^2}$$
Boundary Conditions:
- $T(0) = T_\infty$ ($x=0$)
- $\dfrac{dT}{dx}\big|_{x=L/2} = 0$ ($x=\frac{L}{2}$)
Notes:
- We are going to solve the boundary value problem for $0\leq x \leq \frac{L}{2}$.
- The solution is obtained with the help from a solver BVP from the python library
scipy.integrate.solve_bvp
1. Import necessary libraries¶
Firstly we need to import necessary libraries to this notebook.
numpyis needed for math operations such as numeric arrays, vectorized evaluation of the solution and parameters.scipyfor the core boundary value problem solver bvp.matplotlibfor plotting the results.
import numpy as np
from scipy.integrate import solve_bvp
import matplotlib.pyplot as plt
2. Physical Parameters¶
Below are the physical constants and parameters that define the problem.
These describe the thermal, electrical and geometric characteristics of the current-carrying wire.
# -----------------------------
# Thermal and electrical properties
k = 72.0 # Thermal conductivity [W/(mK)]
h = 2000.0 # Convective heat transfer coefficient [W/(m^2K)]
eps = 0.1 # Surface emissivity
sigma = 5.67e-8 # Stefan–Boltzmann constant [W/(m^2K^4)]
rho_e = 32e-8 # Electrical resistivity [omega m]
# -----------------------------
# Electrical and ambient conditions
I = 2.0 # Electric current [A]
T_inf = 300.0 # Ambient temperature [K]
# -----------------------------
# Geometric parameters
D = 7.62e-5 # Wire diameter [m]
L = 4.0e-3 # Total wire length [m]
A = np.pi * D**2 / 4.0 # [m^2]
3. Auxiliary Constants¶
To simplify the differential equation, we define some auxiliary constants below. These parameters remain constant trhoughout the domain and allow us to write the differential equation in a more compact form.
C1 = 4*h / (k*D)
C2 = 4*eps*sigma / (k*D)
C3 = I**2 * rho_e / (k*A**2)
print(f"C1 = {C1:.3e}, C2 = {C2:.3e}, C3 = {C3:.3e}")
4. Define the Differential Equation and Boundary Conditions¶
Second-Order ODE to First-Order System¶
We started from a second-order differential equation describing the steady-state temperature variation along the wire. $$ \frac{d^2T}{dx^2} - \frac{4h}{kD}(T - T_\infty) - \frac{4\varepsilon\sigma_{SB}}{kD}(T^4 - T_\infty^4) = -\frac{I^2 \rho_e}{k\left(\frac{\pi D^2}{4}\right)^2}. $$
To make it more suitable for numerical solution using solve_bvp, we converted this equation into a system of first-order equations by introducing two new variables:
$$ y_1(x) = T(x), \quad y_2(x) = \frac{dT}{dx}. $$
Substituting these definitions gives:
$$ \begin{cases} y_1' = y_2, \\[4pt] y_2' = C_1 (y_1 - T_\infty) + C_2 (y_1^4 - T_\infty^4) - C_3, \end{cases} $$
where the constants (C_1), (C_2), and (C_3) are defined as above.
Boundary Conditions¶
Because the temperature distribution is symmetric about the midpoint of the wire, we only need to solve for half the wire length, from (x=0) to (x=L/2).
$$ \begin{cases} y_1(0) = T_\infty = 300\text{ K}, \\[4pt] y_2(L/2) = 0. \end{cases} $$
- The Dirichlet boundary condition fixes the temperature at the wire surface ($x=0$).
- The Neumann boundary condition imposes zero temperature gradient at the midpoint ($x=\frac{L}{2}$), representing symmetry and no heat flux through that plane.
Figure 2 – EDO Resolution
#Function defining the ODE system
# y[0] = T
# y[1] = dT/dx
def fun(x, y):
T = y[0]
dTdx = y[1]
d2Tdx2 = C1*(T - T_inf) + C2*(T**4 - T_inf**4) - C3
return np.vstack((dTdx, d2Tdx2))
# Boundary conditions
def bc(ya, yb):
# ya -> values at x=0; yb -> values at x=L/2
return np.array([
ya[0] - T_inf, # T(0) = T_inf
yb[1] # dT/dx(L/2) = 0
])
5. Initial Guess¶
The boundary value solver requires an initial estimate for both the temperature $T(x)$ and its derivative $\frac{dT}{dx}$ along the domain.
As a physically reasonable starting point, we assume:
- The wire is initially at the ambient temperature $(T_\infty = 300\,\mathrm{K})$.
- There is no temperature gradient (flat profile).
This provides a simple and stable initial condition for the iterative solver: $$ T(x) = T_\infty, \quad \frac{dT}{dx} = 0. $$
#Initial guess
x = np.linspace(0, L/2, 50)
y_init = np.vstack((
T_inf * np.ones_like(x), # palpite inicial T(x)=T_inf
np.zeros_like(x) # dT/dx = 0
))
6. BVP Solver¶
Similar to the MATLAB bvp4c function, we use the solve_bvp function from the scipy.integrate library to solve the boundary value problem numerically.
sol = solve_bvp(fun, bc, x, y_init, tol=1e-6, max_nodes=10000)
if sol.success:
print("Solution converged.")
else:
print("Solution did not converge.")
7. Plot the results¶
Finally, we plot the temperature distribution along the wire length to visualize the results.
x_plot = np.linspace(0, L/2, 200)
T_plot = sol.sol(x_plot)[0]
# Max temperature at the center (x=L/2)
T_max = T_plot[-1]
# Tentar usar Plotly para interatividade; se não estiver disponível, usar Matplotlib como fallback
try:
import plotly.graph_objects as go
hover = 'x: %{x:.2f} mm<br>T: %{y:.2f} K'
fig = go.Figure()
fig.add_trace(go.Scatter(
x=x_plot*1e3,
y=T_plot,
mode='lines',
name='T(x)',
line=dict(color='firebrick', width=3),
hovertemplate=hover
))
fig.add_trace(go.Scatter(
x=[x_plot[-1]*1e3],
y=[T_max],
mode='markers+text',
name='T_max',
marker=dict(size=8, color='royalblue'),
text=[f'T_max = {T_max:.2f} K'],
textposition='top right',
hovertemplate='x: %{x:.2f} mm<br>T_max: %{y:.2f} K'
))
fig.update_layout(
title='Distribuição de Temperatura ao longo do fio (0 ≤ x ≤ L/2)',
xaxis_title='x [mm]',
yaxis_title='Temperatura T [K]',
template='plotly_white',
width=800,
height=500
)
fig.update_xaxes(showgrid=True)
fig.update_yaxes(showgrid=True)
fig.show()
except Exception as e:
# Fallback para Matplotlib
plt.figure(figsize=(8,5))
plt.plot(x_plot*1e3, T_plot, 'r', lw=2)
plt.xlabel('x [mm]')
plt.ylabel('Temperatura T [K]')
plt.title('Distribuição de Temperatura ao longo do fio (0 ≤ x ≤ L/2)')
plt.grid(True)
# Anotar a temperatura máxima
plt.scatter([x_plot[-1]*1e3], [T_max], color='blue')
plt.annotate(f'T_max = {T_max:.2f} K', xy=(x_plot[-1]*1e3, T_max),
xytext=(x_plot[-1]*1e3*0.9, T_max + 2),
arrowprops=dict(arrowstyle='->', color='black'))
plt.show()
print(f"T_max no centro (x=L/2): {T_max:.2f} K")
8. Conclusion¶
The steady‑state axial temperature distribution of the current‑carrying bare wire was obtained by solving the boundary value problem on half the length $(0 \leq x \leq \frac{L}{2})$ using symmetry $(T(0)=T_\infty$ and $\frac{dT}{dx}(L/2)=0)$.
The solution rises monotonically from the cooled end to a maximum at the center, where the axial heat flux vanishes. For the given parameters the computed centerline maximum temperature is $T_{max} = 781.60 \, K$.
Reformulating the second‑order ODE as a first‑order system enabled robust convergence with solve_bvp, yielding a smooth temperature field suitable for further parametric or interactive exploration.